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Solving Navier-Stokes Equations with 
Advanced Turbulence Models 
on Three-Dimensional Unstructured Grids 

Qunzhen Wang* Steven J. Massey^ and Khaled S. Abdol-Hamid* 
Analytical Services & Materials, Inc., Lancaster, CA 
Neal T. Frink fj 

NASA Langley Research Center, Hampton, VA 

USM3D is a widely- used unstructured flow solver for simulating inviscid and viscous 
flows over complex geometries. The current version (version 5.0) of USM3D, however, 
does not have advanced turbulence models to accurately simulate complicated flows. 

We have implemented two modified versions of the original Jones and Launder k-e two- 
equation turbulence model and the Girimaji algebraic Reynolds stress model in USM3D. 

Tests have been conducted for two flat plate boundary layer cases, a RAE2822 airfoil 
and an ONERA M6 wing. The results are compared with those of empirical formulae, 
theoretical results and the existing Spalart-Allmaras one-equation model. 


Introduction 

T HE unstructured-grid methodology offers some 
significant advantages compared to the tradi- 
tional structured-grid method for simulating flows over 
complex aerodynamics shapes. This is mainly due 
to the promise that the construction of an unstruc- 
tured grid around complex configurations such as an 
aircraft requires much less time than a comparable 
block-structured grid. Furthermore, local refinement 
of unstructured grids can be carried out more eas- 
ily to improve the accuracy of the simulation. While 
more work remains to be done to fully realize their 
potential, much progress has been made in modeling 
complicated flows on unstructured grid (see Mavriplis 1 
for a review) . 

One important phenomena for complex viscous flows 
is turbulence, which is difficult to simulate due to the 
existence of a wide range of scales. There are many 
types of method to deal with turbulence, ranging from 
the simplest algebraic model to the most accurate di- 
rect numerical simulation. For most of the turbulence 
models, the Reynolds stress is assumed to be related 
to the mean strain rate by the eddy viscosity. Such 
turbulence models may be divided into zero-equation 
model (i.e., algebraic model), one-equation model, and 
two-equation model depending on the number of trans- 
port equations need to be solved to obtain the eddy 
viscosity. The Reynolds stress model does not use the 
concept of eddy viscosity. Instead, a transport equa- 
tion for each component of the Reynolds stress tensor 
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is solved directly. While the above models only solve 
the mean flow, large eddy simulation (LES) solve the 
large scale fluctuations in addition to the mean flow 
and only the effect of small scales (i.e., subgrid scales) 
are modeled. Finally, the most accurate method is 
direct numerical simulation (DNS), where both mean 
flow and all the fluctuations are solved directly. 

While the Reynolds stress model, LES and DNS 
methods are much more accurate than the eddy vis- 
cosity based method, they require prohibitive amounts 
of CPU time and memory. Therefore, the most 
widely used turbulence models in industry are still 
those based on the concept of eddy viscosity. This 
is especially true for unstructured grid CFD codes. 
For example, the predominantly utilized turbulence 
model in the finite-element unstructured grid code of 
Mavriplis 2 is the algebraic model of Baldwin and Lo- 
max 3 although extension has been made to include a 
two-equation model (Mavriplis and Martinelli 4 ). The 
unstructured nodal-based finite volume code FUN3D 
(see Anderson, 5 Anderson and Bonhaus 6 ), contains 
two one-equation models, one by Baldwin and Barth' 
and the other by Spalart and Allmaras. 8 

USM3D is a tetrahedral cell-centered unstructured 
flow solver for simulating inviscid and viscous flows 
over complex geometries. It is developed by Frink 9-12 
at NASA Langley Research Center and is now be- 
ing widely used in both industry and government. 
This code is part of the TetrUSS flow analysis and 
design package which won the 1996 NASA Software 
of the Year award. The algorithm of USM3D con- 
sists of a cell-centered, upwind-biased, finite-volume, 
implicit/explicit algorithm capable of solving the com- 
pressible Euler and Navier-Stokes equations on un- 
structured tetrahedral meshes. Like most other un- 
structured CFD codes, however, USM3D does not 
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have advanced turbulence models to accurately simu- 
late complex flows. The current production version of 
the code only has the one-equation turbulence model 
developed by Spalart and Allmaras, 8 although a two 
equation k-e model was implemented in a previous re- 
search version of USM3D (Kwon and Hah 13 ). It is well 
known that one-equation turbulence models are not 
adequate for complex flows such as separated flows or 
shear flows. The objective of this study is to compare 
the performance of the recently added advanced tur- 
bulence models (see Wang et al . 14 ) with experimental 
data, theoretical results and the results from the exist- 
ing Spalart- Allmaras one-equation model implemented 
by Frink. 12 


Governing Equations 

The integral form of the Navier-Stokes equations, 
which govern the compressible, Newtonian flow of fluid 
in the absence of external forces, can be written as 


^ J JJ CldV + I J F ^' fldS= j j G (Q )-nd,S 

(1) 

where the first term (time changing rate) is integrated 
over the volume of a bounded domain while the sec- 
ond term (convection) and the third term (diffusion) 
are integrated over the boundary of this domain. The 
quantities Q, F(Q), and G(Q) are vectors with five 
components and Eq. (1) contains five equations corre- 
sponding to the conservation of mass, momentum, and 
energy. The unknowns in Eq. (1) are 


Q 


p 

pu 

pv 

pw 

e 


(2) 


where p , u, v, w and e are density, three velocity com- 
ponents, and total energy, respectively. The inviscicl 
flux is 


F(Q) • fi = 
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pv 
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h z 
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( 3 ) 


where p is the pressure and h x , h y , and h z are Carte- 
sian components of the exterior surface unit normal h 
on the boundary of the domain. The viscous flux is 


G(Q) • h 


Mo o 

Re L 


{h x G i + n y G2 + fi z Gz) 


( 4 ) 


where 


Gi = 


0 

Til 

Ta 

Ta 

'U'jl’ij Qi 


( 5 ) 


Moo is the free-stream Mach number and ReL is the 
Reynolds number based on a typical length scale (e.g., 
the total length of the flat plate) . 

The total stress T t] and heat flux qi can be divided 
into a laminar part (denoted by superscript L) and a 
turbulence part (denoted by superscript T) 

■/ = "/; + t'i (6) 

Qi = Qi + Qi ( 7 ) 

with the laminar part being 



2 du k 
3' ij dx k 


(8) 


Qi = - 


p L dT 
(7 - l)Pr L dxi 


(9) 


where p L is the molecular viscosity, T is the tem- 
perature, Pr L is the molecular Prandtl number, and 
7 = 1.4 is the ratio of specific heats. Following Eqs. (8) 
and (9), the Reynolds stress and turbulent heat flux 
can be approximated as 



2 du k 

3 ij dx k 



(10) 


p T dT 
(7 - l)Pr T dxi 


( 11 ) 


Note that the last term in Eq. (10) is presented only 
when a transport equation for turbulent kinetic energy 
k is solved. 


Transport Equations for k and e 

In the k-e model, two transport equations, one for 
turbulent kinetic energy k and the other for dissipation 
rate e, are solved and the eddy viscosity p T in Eqs. (10) 
and (11) are then calculated based on k and e 


p T = CpUptj (12) 

where is a damping function and Gy = 0.09. The 
transport equation for k is 


dpk dpkuj d 
~dt ^ 


where 


dxj dxj 


Pk 


dk 

dxj 
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S k = P^-p(l+T)e 
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ReL 


= S k 
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P = rl 


duj 
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j P 

Pk = P H cr fe = 1.0 

& k 


(13) 

(14) 

(15) 

(16) 


Similarly, the transport equation for e can be written 
as 


dps 

~dt 


dpsuj 

dxj 


d 

dxj 


ds ] Mpo _ „ 

^ e dxj ReL 


(17) 
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L k = 2 p L 
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(19) 

T 

p e = p L +^— cr e = 1.3 <761 = 1.44 C e2 = 1.92 

c e 

(20) 

ok 2 

h = 1 - 0-3 exp(— i? 2 ) Rt = ^j- (21) 
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In all of the equations, as well as, all the figures shown 
in this paper, unless explicitly stated otherwise, the 
length is normalized by a characteristic length L, the 
velocity by a the density by poo, the viscosity by 
Poo, the turbulent kinetic energy k by o^,, and the 
dissipation rate £ by poo ato / where doc is the free- 
stream speed of sound, poo is the free-stream density, 
and poo is the free-stream molecular viscosity. 

The compressibility correction T in Eq. (14) and 
the damping function in Eq. (12) can take differ- 
ent forms. The two most widely used compressibility 
corrections are Sarkar et al. 15 model 


r = m 2 


(22) 


and Wilcox 16 model 


T = (M 2 - M t 2 0 ) H(M t - M t0 ) (23) 


where H{ x) is the Heaviside step function, the turbu- 
lent Mach number M t = \fk/a with a being the local 
speed of sound. The damping function could take one 
of the following three forms: a) Jones and Launder 17 
form 


fn = exp 


3.41 

(i + S ) 2 


(24) 


b) Van Driest form (Nagano and Hishida 18 ) 


fn = 1 - exp 



(25) 


c) Speziale et al. 19 ) form 


fn = 



tanh 



(26) 


For the results shown in this paper, ivisc=6 refers to 
the Jones and Launder form of the damping function, 
Eq. (24), with no compressibility correction (i.e., T = 
0) and ivisc=7 refers to the Jones and Launder model 
modified by Carlson 20 as given in Eqs. (27-31). 


= 0.41 a = 1.15 C e 2 = 1.9 


(27) 


C el = a [ 1 + 0.2174- 


(28) 


1 5\ 

1 + 0.2174— 


(C-62-Cry^ 



fn = exp 


6 


(29) 

(30) 

(31) 


Algebraic Reynolds Stress Model 

The k-e model discussed above is the standard 
model, which is also called the “linear model” be- 
cause the turbulent stress is linearly related to the 
mean strain rate by the eddy viscosity as is clear from 
Eq.(10). However, various direct numerical simulation 
(DNS) data have shown that the turbulent stress does 
not vary linearly with the mean strain rate. A more 
accurate model is the Reynolds stress model where a 
transport equation for each component of the Reynolds 
stress tensor is solved directly and the concept of eddy 
viscosity is not used. However, the Reynolds stress 
model requires a tremendous amount of CPU time and 
memory and, thus, is seldom used in real complex en- 
gineering applications. Something in between is the 
algebraic stress model, in which a nonlinear term is 
added to the turbulent stress (so it is also called the 
“nonlinear model”). Following Girimaji, 21 the turbu- 
lent stress in the algebraic stress model is 


T T 

T H = V 


Sin ~ ~Si 


1 „ dui 


dx k 


2 pkSij 


+2p T K 1 - [S lk W kj - W lk S kj ] 


+2p T I< 1 - 


Si k S k j S k [S[ k 


(32) 


where the mean strain rate and mean vorticity are 

(33) 


_ 1 ( dui t duj 


dxj dxi 


rrr 1 ( 9Ui dUj 

Wii = ~ 1 


2 \ dx j dxi 
respectively. In Eq.(32), I\\ and K 2 are given by 


(34) 
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G 2 
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and G 2 and G 3 are 
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(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 

(52) 

(53) 

(54) 

(55) 


Furthermore, instead of G p = 0.09 for the linear 
model, G m = — G\ is applied in the nonlinear model. 


Numerical Procedure 

The details of the numerical procedure for solv- 
ing equation (1) are discussed in Frink 10-12 and only 
a brief overview is given here. The spatial compu- 
tational domain is divided into a finite number of 
tetrahedral cells and a finite-volume discretization is 
applied to each cell. This procedure results in a set 
of volume-averaged state variables Q which are in bal- 
ance with the area-averaged inviscicl flux F(Q) and 
viscous flux G(Q). Inviscicl fluxes are obtained across 
each cell face using either the Roe 22 flux-difference 
splitting approach or the Van Leer 23 flux-vector split- 
ting technique. The data at nodal points could be ob- 
tained from the cell-averaged data by either an inverse- 
distance weighted averaging scheme or a Laplacian- 
weighted averaging scheme. The viscous fluxes are 
approximated at the cell-face centroids by linear re- 
construction. An implicit time integration algorithm 
using the linearized backward Euler time differencing 
approach is applied to update the solution. The result- 
ing linear system of equations are solved at each time 
step with a subiterative procedure by a point- Jacobi 
method. Convergence to the steady state solution is 
accelerated by advancing the equations at each cell in 
time by the maximum permissible time step in that 
cell. The CFL number is scaled according to the de- 
viation of cell aspect ratio from the ideal value of an 
isotropic tetrahedron. 

The solution procedure for the two turbulence trans- 
port equations (13) and (17) are similar to that for 
solving the Navier-Stokes equations. Equations (13) 
and (17) are solved separately from the flow govern- 
ing equations and from each other using the same 
backward Euler time integration scheme. The k and 
e equations can be solved using either first-order or 
second-order schemes. For the second-order method, 
either Roe’s SuperBee limiter or the Minimum Mod- 
ulus (Min-Mod) limiter could be applied. To allow a 
large CFL number, implicit method is used to solve 
the k and £ equations. 

Two input parameters ko and /ij are needed to spec- 
ify the initial conditions and limit the smallest values 
of k and e. The initial conditions are k = ko and 
£ = £0 = G^pfcp/ /j,Q . The turbulent kinetic energy and 
dissipation rate are not allowed to be smaller than ko 
and £ 0 , respectively. On solid surfaces, the boundary 
conditions are k = ko and £ = / p (M^ / Re^) 2 . Far 

field boundary conditions are applied by extrapolating 
k and e from the interior for outflow boundaries and 
taken from the free-stream for the inflow boundaries. 

The k-e model has been coupled with wall function 
formulations to reduce the need for resolving the flow 
in the near wall region. The following three differ- 
ent wall functions could be used together with the k-e 
model: (1) iwallf=l: the original wall function devel- 
oped by Frink ; 12 (2) iwallf=2: a wall function similar 
to that used in PAB3D (see Abdol-Hamid et al. 2A for 
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details); (3) iwallf=3: same as the original wall func- 
tion except the velocity components are zero at the 
solid surface. For iwallf=l, a slip velocity boundary 
condition is obtained by solving the Spalding formula 


y + = u + +e~ KB 


:x+ + (kw +) 2 

- 1 - ki C - - — — — 


( ku + ) 3 


6 

(56) 

with k = 0.4 and B = 5.5 using Newton-Raphson iter- 
ation method while the velocity at the solid boundary 
for iwallf=2 and iwallf=3 are zero. 


Results 

Flat Plate Cases: BLT2 and BLT3 

Two tetrahedral grids for a simple flat plate bound- 
ary layer have been generated using the grid generator 
VGRID (Pirzadeh 25 ). For both grids, results of the 
newly added k — e models are compared with the ex- 
isting Spalart-Allmaras model and with experimental 
and theoretical data. The computational domain is 
from x = —0.5 to x = 1.0 in the streamwise direction, 
from y = —0.02 to y = 0.02 in the spanwise direc- 
tion (—0.05 < y < 0.05 for BLT2), and from z = 0 to 
z = 0.22 in the wall normal direction. The free-stream 
Mach number and Reynolds number are = 0.5, 
ReL = 2 x 10 5 6 , respectively. The grid characteristics 
for the flat plate cases are summarized in Table 1. 


Table 1 Characteristics of flat plate grids. 


Case 

Cells 

Nodes 

Faces 

B.Nds 

B.Fcs 

BLT2 

BLT3 

48,497 

37,483 

9,629 

8,038 

99,805 

78,328 

2,813 

3,364 

5,612 

6,724 


The grids for the flat plate cases are shown in Fig- 
ure 1 . There are four cells in the spanwise direction for 
BLT2 and two cells for the BLT3 grid. Near the wall, 
the grid spacings for the BLT3 grid are much smaller 
than that for BLT2 but the grid spacings for BLT3 is 
much larger than that for BLT2 far away from the wall. 
For the BLT3 grid, the first node point away from the 
wall has a y + ss 1.8, while for the BLT2 grid y + « 223. 
Therefore, a wall function has to be applied for BLT2 
while no wall function is needed for the BLT3 grid. 

BLT2 Flat Plate - Wall Function Case 

Results using the Spalart-Allmaras and linear k-e 
models are shown in Figure 2. The residual history for 
the Spalart-Allmaras and Carlson modified k-e model 
both converge rapidly at nearly the same rate to a 
level of 9 orders of magnitude smaller than the initial 
residual. The first linear model (ivisc=6) converged 
much slower after the 500 </l time step and only de- 
creased by less than 4 orders of magnitude after 4,000 
time steps. The CFL number is allowed to increases 
dynamically from 1 to 200 according to the residual: it 
increases when the residual is decreasing and decreases 
when the residual is increasing. For each of the three 



a) Full BLT2 grid. 


b) Wall of BLT2 grid. 



c) Full BLT3 grid. 


d) Wall of BLT3 grid. 


Fig. 1 Whole domain and near wall regions of BLT2 
and BLT3 grids for a flat plate boundary layer. 


models tested, the CFL number rapidly increased to 
its maximum value within 40 time steps and remained 
at 200 for the duration of the run, indicating a robust 
convergence. The velocity profile is compared with 
the empirical formula of Spalding, Eq. (56), while the 
skin friction coefficient is compared with the theoreti- 
cal values for fully turbulent flow. 

C f = 0.0583(Re x )“ 1/5 (57) 


It is evident that u + from the second (ivisc=7) k- 
e model is the closest to Spalding’s formula with the 
Spalart-Allmaras results the next closest. The first 
k-e model initially matches the data as dictated by 
the wall function, but then significantly under pre- 
dicts u + . The skin friction coefficient from the first 
k-e model is significantly larger than the theoretical 
value whereas the Spalart-Allmaras model gives a re- 
sult slightly smaller than the theoretical value and the 
second k-e model initially matches theory closely be- 
fore predicting a slightly larger value. In summary, 
for the coarse flat plate grid using a wall function, 
the Carlson modified linear k-e model performed best, 
with the Spalart-Allmaras model nearly as good and 
the first k-e model not in good agreement with empir- 
ical or theoretical data. The CPU time per timestep 
per cell on an Intel Pentium II 300 MHz was 229 /zsec 
for the Spalart-Allmaras model and 13% more for the 
k-e models. 


BLT3 Flat Plate - Grid Resolved Case 

The results from the Spalart-Allmaras and k-e mod- 
els are shown in Figure 3. The residual history for all 
three models converge rapidly at roughly the same rate 
to a level of 5 orders of magnitude smaller than the 
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Log Residual for BLT2 Case 




a ) 


b) 


u + vs. y + BLT2 Case 



c ) 


Skin Friction BLT2 Case 



d) 


Fig. 2 History of a) Log residual and b) CFL num- 
ber. Converged c) streamwise velocity component 
in wall coordinates at x/c = 0.5 and d) skin friction 
coefficient. BLT2 case. 


initial residual. The temporary hump in the residual 
trace around n = 1,000 is characteristic of the estab- 
lishment of turbulence for a non- wall function grid. 
For each of the three models tested, the CFL number 
rapidly increased to its maximum value within 75 time 
steps and remained at 200 for the duration of the run, 
indicating a robust convergence. 

Inspection of the velocity profile and skin friction 
plots, see Figure 3, confirms that for the BLT3 grid a 
wall function is not needed. It is evident that u + from 
the Spalart-Allmaras model is closer to the Spalding 
curve with the second k-E model under predicting u + 
slightly and the first k-E model under prediction sig- 
nificantly. The skin friction coefficient from the first 
k-E model is significantly larger than the theoretical 
value whereas the Spalart-Allmaras model gives a re- 
sult slightly smaller than the theoretical value and the 
second k-E model matches theory the closest predicting 
a slightly larger value. 

In summary, for the BLT3 refined flat plate grid 
without using a wall function, the Carlson modified 
k-E model performed best for predicting skin friction 
and is a close second to Spalart-Allmaras for veloc- 
ity profile. Experience with the linear models in the 
structured code PAB3D 24 indicates that with grid re- 
finement the accuracy of the linear models will improve 
significantly. The CPU time per timestep per cell 
on an SGI with two IP30 processors was 194 /nsec for 
the Spalart-Allmaras model and 10% more for the k-E 
models. 


Log Residual for BLT3 Case 



a ) 



b) 


u + vs. y + BLT3 Case Skin Friction BLT3 Case 



Fig. 3 History of a) Log residual and b) CFL num- 
ber. Converged c) streamwise velocity component 
in wall coordinates at x/c = 0.5 and d) skin friction 
coefficient. BLT3 case. 

RAE2822 Airfoil 

In this section results for the previously mentioned 
turbulence models of Spalart-Allmaras model, k-E 
(ivisc=6), k-E (ivisc=7) will be compared with the 
Girimaji algebraic Reynolds stress model (ARSM) as 
well as the experimental results of case 10 of Cook et 
al. 26 for a transonic airfoil. 

The computational domain extends 6 chord lengths 
away from the airfoil in all directions with a grid width 
0.2c spanning 2 cells. The free-stream Mach number 
of the boundary layer, Reynolds number and corrected 
free air angle-of-attack is M 0 c = 0.75, Re = 6.2 x 10 6 
and a = 2.81° respectively. The grid consist of 29,976 
cells, 8,477 nodes, 66,772 faces, 6,820 boundary nodes 
and 13,640 boundary faces. At the first node away 
from the wall y + ~ 0.8. Full and closeup views of the 
grid are shown in Figure 4. 

Unlike the previous flat plate cases, the residual and 
CFL history (Figure 5) are much more noisy, which 
may be due to unsteadiness in the region behind the 
shock. All of the models converge between 3 and 3.5 
orders of magnitude after 4,000 time steps. Initially 
the ARSM shows a great deal of oscillation since it 
is more dependent on the initial condition. Also ob- 
served is a large increase in residual for the second 
linear model before finally settling down. The CFL 
history reveals that for the first 400 time steps each 
model is behaving similarly with the exception of the 
ARSM whose CFL history appears to be shifted to 
the left. After n = 400 the CFL number for all models 
oscillate periodically at high frequency. 

From the lift and drag coefficient history plots (Fig- 
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ure 5) it is evident that the Spalart-Allmaras, first 
k-e and Girimaji models follow the same convergence 
trend while the second linear model oscillates much 
more before settling down at around n = 2, 500. The 
results at n = 4, 000 along with experimental data 
and results from PAB3D for the first linear model are 
compared in Table 2. Final results show that the 
Spalart-Allmaras model is in closest agreement with 
experimental data. The next best results for USM3D 
are produced by the first linear model, the ARS model 
and finally the second linear model. As an indicator 
of grid dependence for the advanced models, results in 
good agreement with experimental data from PAB3D 
are shown for the first linear model. 


Table 2 Comparison of Cl and Cd with experi- 
mental data of Cook et al. and PAB3D linear k-e 
model (ivisc=6). 24 


Case 

C L 

% diff 

C D 

% diff 

Experiment 

0.743 

- 

0.0242 

- 

Spalart-Allmaras 

0.732 

-1.5 

0.0214 

-11.6 

k-e (ivisc=6) 

0.691 

-7.0 

0.0199 

-17.8 

k-e (ivisc=7) 

0.649 

-12.7 

0.0178 

-26.4 

Girimaji ARSM 

0.803 

8.1 

0.0293 

21.1 

PAB3D k-e 

0.720 

-3.1 

0.0257 

6.2 


In Figure 6 the predicted coefficient of pressure from 
each model is plotted with the experimental data of 
Cook. 26 All models miss the location of the first suc- 
tion peak by approximately 2%c. This is likely due to 
the fact that the current implementation of the turbu- 
lence models in USM3D does not allow for setting a 
trip location for the boundary layer, while in the exper- 
iment it was set to 3%c. Other areas of disagreement 
with experiment and among the models themselves is 
the shock location and region aft of the shock. The 
prediction of shock location by the second linear model 
was in excellent agreement while the other models all 
predicted different locations-all aft of the experimen- 
tal location. Following the shock, all models predicted 
a lower Cp than experiment, with the Girimaji model 
in the best agreement. 

The CPU time per timestep per cell on an Intel Pen- 
tium II 300 MHz was 207 p,sec for the Spalart-Allmaras 
model, 224 y^sec for the first k-e model, 234 /xsec for the 
second k-e model and 302 /rsec for the Girimaji model. 

ONERA M6 Wing 

The tetrahedral viscous grid for the ONERA M6 
wing was generated using VGRID and is similar in 
construction to those in Frink. 12 The grid consist 
of 338,417 cells, 59,496 nodes, 682,257 faces, 5,425 
boundary nodes and 10,846 boundary faces. On the 
wing surface y + « 2 for M x = 0.8447, Re mac = 
11.7 x 10 6 and a = 5.06°. The computational do- 
main is bounded by a rectangular box defined by 



Fig. 4 RAE2822 airfoil grid, a) whole domain; b) 

near airfoil. 


Log Residual for RAE2822 



CFL History for RAE2822 

300 

A Spalart-Allmaras 

250 - B Linear k-e (ivisc=6) 

0 Linear k-e (ivisc=7) 



n 


b) 



C D History for RAE2822 


Experiment 
Spalart-Allmaras 
Linear k-e (ivisc=6) 
Linear k-e (ivisc=7) 
Girimaji ARSM 



c) d) 

Fig. 5 History of a) Log residual, b) CFL number, 
c) Cl and d) Cd for RAE2822 airfoil. 


Cp RAE2822 Airfoil 



Fig. 6 Comparison of coefficient of pressure from 
each model (n = 4, 000) with the experimental data 
of Cook 26 . 
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—6.5 < x < 6.5, 0 < y < 4, and —6.5 < z < 6.5, in 
aerodynamic coordinates relative to a semispan length 
of 1. Surface and symmetry plane meshes are shown 
in Figure 7. 

For this ONERA M6 wing case, results are presented 
for the Spalart-Allmaras model and the second lin- 
ear (ivisc=7) model with varying numbers of Jacobi 
iterations in the solution of the k-e equations. In Fig- 
ure 8, the residual history is seen to be fairly similar 
among the models with the Spalart-Allmaras being the 
smoothest and both linear cases being more noisy un- 
til n = 500. To aid solution stability, the original CFL 
ramping scheme in USM3D was used to bring the CFL 
number to a constant value of 75, rather than letting 
it vary dynamically up to 200 as in previous cases. 

The coefficient of lift and drag histories (Figure 8) 
confirm that the reduction in the number of iterations 
on the k-e equations causes no inaccuracy in the solu- 
tion. This is significant because the reduction brings 
a 10% time savings. While the two linear cases agree 
with each other, the Spalart-Allmaras model predicts 
a 4% lower Cl while agreeing exactly with the linear 
models prediction of Cd ■ The downward trend in both 
Cl and Cd plots indicate a slight lack of convergence 
despite the drop in log residual of over 3 orders of mag- 
nitude. This may be a result of the unsteadiness of the 
flow field. 

Limiting surface streamlines simulating “oil-flow” 
patterns shown in Figure 9 show a significant amount 
of shock separated flow beyond the 77 ~ 0.65 for 
both models, with the Spalart-Allmaras model results 
showing the strongest separation. The Cp plots in 
Figure 10 demonstrate this with the exception of the 
77 = 0.90 plane where Spalart-Allmaras matches the 
shock location well. It should be noted however, that 
the solutions in the tip region, 77 > 0.90, can be partic- 
ularly sensitive to a variety of factors such as grid den- 
sity and turbulence model (see Rumsey and Vatsa 27 ), 
and should be explored further in future studies. 

The CPU time per timestep per cell on an SGI 
with two IP30 processors was 216 /xsec for the Spalart- 
Allmaras model, 223 /i sec for the second (ivisc=7) k-e 
model with 2 Jacobi iterations (nstagek=2) on the k- 
s equations and 249 fj,sec for the second linear model 
with 6 Jacobi iterations. 

In summary, the k-£ model is seen to perform well in 
a complex three-dimensional flow with only a small in- 
crease in CPU time required over the existing Spalart- 
Allmaras model. 

Concluding Remarks 

A systematic study has been conducted to as- 
sess the accuracy of two newly implemented turbu- 
lence models; modified versions of the standard linear 
two-equation k-e model and the non-linear algebraic 
Reynolds stress model of Girimaji. Initial results of 
tests on the flat plate, airfoil, and wing cases indicate 



a ) b) 


Fig. 7 Surface mesh for a) entire wing and b) root 
region of symmetry plane. 


Log Residual for ONERA M6 CFL History OM6 Wing 




C L History OM6 Wing C D History OM6 Wing 




Fig. 8 History of a) Log residual, b) CFL number, 
c) Cl and d) C D for ONERA M6 wing. 


OM6 USM3Dns Solution - 2000 cycles 

M=0.8447, AOA=5.06. Re=11.7M 



Fig. 9 Comparison of k-e (ivisc=7) and Spalart- 
Allmaras model computed surface “oil-flow” pat- 
terns and flooded pressure contours for the ON- 
ERA M6 wing. 
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Cp 0M6 Wing r|=0.44 C p 0M6 Wing r|=0.65 




C p 0M6 Wing r|=0.80 C p 0M6 Wing r|=0.90 




Fig. 10 Comparison of coefficient of pressure from 
each model (n = 2, 000) with the experimental data 
of Schmitt 28 . 

that the new two-equation models yield comparable 
accuracy and efficiency to that of the Spalart-Allmaras 
one-equation model. Work is currently underway to 
further examine factors such as grid sensitivities and 
solution details. 
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